Spatial and Temporal Variations in Microbial Diversity on Foxfonna
Archana Dayal & André Soares
11 May, 2021
- Loading libraries
- Importing data
- Testing batch effects
- 1 - Library sizes
- 2 - Evaluating dataset contamination with
decontam - 3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
- 3b - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Snowpack_Layer
- 3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month
- 4 - Top SVs and their sequences
- 5a - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (no filter)
- 5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)
- 5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)
- 6 - Unconstrained ordination
- 7 - Alpha-diversity metrics.
- All measures for Sample_ID, no stats, facet by DNA/cDNA
- All measures by Sample_ID
- All measures for Month, stats included, facet by DNA/cDNA
- All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
- All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
- All measures for Fox_Aspect, no stats, facet by DNA/cDNA
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed
- All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed
- All measures for Snowpack_Layer, no stats, facet by DNA/cDNA
- All measures for Snowpack_Layer, stats included, facet by DNA/cDNA
- All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA
- All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed
- All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed
- All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed
- 8 - Alpha diversity with rarefaction to the smallest number of reads in any sample
Loading libraries
library(dada2); packageVersion("dada2")## [1] '1.18.0'library(ggplot2); packageVersion("ggplot2")## [1] '3.3.3'library(phyloseq); packageVersion("phyloseq")## [1] '1.34.0'library(data.table); packageVersion("data.table")## [1] '1.13.6'library(dplyr)
library(plyr)
library(magrittr)
library(RColorBrewer)
library(picante)
library(igraph)
library(cowplot)
library(plotly)
library(lemon)
library(GGally)
library(randomcoloR)
library(knitr)
library(kableExtra)
library(ggpmisc)
library(ggpubr)
library(decontam)
library(stringr)Importing data
SVt = readRDS("../seqtab_final.rds")
# SVt_cnn = collapseNoMismatch(SVt, minOverlap = 20,
# identicalOnly = FALSE, vec = TRUE, verbose = TRUE)
# saveRDS(SVt_cnn, "../seqtab_final_cnn.rds")
# SVt = readRDS("../seqtab_final_cnn.rds")
SVt = SVt[!rownames(SVt) %in% "Undetermined", ]
Tax = readRDS("../tax_final.rds")
metadata<-read.csv("../../220319_for_graph.csv", header=TRUE)
# metadata[, 10:47] <- as.numeric(metadata[, 10:47])
metadata$Snowpack_Layer = gsub("Glacial Ice", "GL ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Sup. Ice", "SUP ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Top", "TOP", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Mid", "MID", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Basal", "BASAL", metadata$Snowpack_Layer)
order<-rownames(SVt)
metadata$Sample_ID=as.character(metadata$Sample_ID)
metadata_ord<-metadata[match(order, metadata$Sample_ID),]
rownames(metadata_ord) <- metadata_ord$Sample_ID
metadata_ord$Sample_Name <- factor(metadata$Sample_Name, levels=unique(metadata$Sample_Name))
ps <- phyloseq(otu_table(SVt, taxa_are_rows=FALSE),
sample_data(metadata_ord),
tax_table(Tax))
ps.b <- prune_taxa(taxa_sums(ps) > 0, ps)
ps.b = subset_taxa(ps.b, Kingdom %in% c("Archaea", "Bacteria"))
ps.b = subset_samples(ps.b, Fox_Aspect != "N_NE_D")Testing batch effects
#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)
readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
samples = rownames(as.data.frame(sample_sums(ps.batch))),
batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]
my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
c("2","3"), c("2","4"), c("2","5"),
c("3","4"), c("3","5"),
c("4","5"))
ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
color=batch, fill=batch)) +
geom_violin(alpha=0.3) +
geom_jitter(shape=16, position = "jitter") +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Batch Number") +
guides(fill=FALSE) +
xlab("Extraction batch") +
ylab("Number of reads obtained") +
theme(axis.text.x = element_text(vjust = 1,
hjust = 1,
size = 12),
strip.text = element_text(size = 12, margin = margin()))
1 - Library sizes
sort(sample_sums(ps.b))## C05-cDNA-29 D12-cDNA-MM C10-cDNA-34 A02-cDNA-MM-H20 C08-cDNA-32
## 0 0 5 7 25
## C07-cDNA-31 B07-cDNA-17 E01-cDNA-MM-H20 B11-cDNA-22 B09-cDNA-19
## 27 43 55 123 187
## D01-cDNA-37 A01-cDNA-MM A11-cDNA-9 D10-cDNA-46 B10-cDNA-21
## 262 290 2320 3844 9170
## A10-cDNA-8 B05-cDNA-15 B08-DNA-4 B11-DNA-7 B08-cDNA-18
## 15595 15781 16781 29350 29808
## C03-cDNA-27 A08-cDNA-6 A06-cDNA-4 A09-cDNA-7 A07-cDNA-5
## 31806 33041 37485 40604 44392
## A11-DNA-9 A05-cDNA-3 A08-DNA-6 B02-cDNA-12 D08-DNA-8
## 50662 52562 53872 57398 57493
## E02-DNA-MM C02-cDNA-26 D07-cDNA-43 D06-DNA-6 D02-DNA-2
## 58537 59570 59805 60653 61230
## A01-DNA-MM B12-DNA-8 B12-cDNA-23 A12-DNA-10 C07-DNA-5
## 61807 62146 62235 63891 64148
## B04-cDNA-14 A07-DNA-5 A02-DNA-MM-H20 A04-cDNA-2 B03-cDNA-13
## 68029 69142 69499 70613 76895
## D09-DNA-9 B02-DNA-12 A05-DNA-3 C04-cDNA-28 C05-DNA-3
## 80028 83134 84732 86983 87568
## C04-DNA-2 D06-cDNA-42 A09-DNA-7 B01-cDNA-11 B04-DNA-14
## 96798 101335 108915 110464 113222
## C03-DNA-B2-1 B03-DNA-13 C08-DNA-6 A12-cDNA-10 B06-cDNA-16
## 114048 118548 122463 123046 123687
## C12-cDNA-36 D02-cDNA-38 D03-DNA-3 C06-DNA-4 E03-DNA-MM-H20
## 123719 137347 138813 139808 140502
## D07-DNA-7 C12-DNA-10 D04-DNA-4 B09-DNA-5 D08-cDNA-44
## 140668 148596 151009 155211 158902
## A06-DNA-4 C01-DNA-9 B10-DNA-6 A10-DNA-8 C09-DNA-7
## 166098 169126 170463 170498 173577
## B01-DNA-11 E01-DNA-4 C06-cDNA-30 C09-cDNA-33 D03-cDNA-39
## 174897 175874 188016 190708 192008
## A03-cDNA-1 A03-DNA-1 B06-DNA-2 C11-cDNA-35 A04-DNA-2
## 193167 196436 198649 198996 200229
## B07-DNA-3 D04-cDNA-40 C10-DNA-8 D11-cDNA-47 D01-DNA-B4-1
## 204577 204668 212618 220186 238199
## D09-cDNA-45 D11-DNA-2 D05-DNA-5 D10-DNA-B5-1 C01-cDNA-25
## 240324 252560 267031 280585 293592
## D05-cDNA-41 D12-DNA-3
## 302761 312011df <- as.data.frame(sample_data(ps.b)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps.b)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_type)) + geom_point()
2 - Evaluating dataset contamination with decontam
#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)
sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
neg="is.neg")
ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")
ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")
ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
ps.final.noncontam.noconts.nodoubs)
ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)
ps.final.r <- transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]ps.final.r.f## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2886 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 2886 taxa by 6 taxonomic ranks ]ps.final## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]sum(rowSums(otu_table(ps.final)))## [1] 6750655gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))Setting order for all plots
sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))3a - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Fox_Aspect, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
ps.final.r.glom <- tax_glom(subset_samples(ps.final.r,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)
ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]
#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)
ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)
names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)
tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")
ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa,
levels = c("Alphaproteobacteria", "Deltaproteobacteria",
"Gammaproteobacteria",
"Other Proteobacteria",
"Actinobacteria",
"Acidobacteria","Bacteroidetes",
"Chloroflexi","Cyanobacteria",
"Firmicutes","Gemmatimonadetes","Other"))ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Month, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill") +
scale_fill_manual(values = tol18rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 12),
strip.background = element_rect(size=1),
legend.text = element_text(size= 12)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance\n")
3b - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Snowpack_Layer
ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Month, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Snowpack_Layer,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill") +
scale_fill_manual(values = tol18rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 11),
strip.background = element_rect(size=1),
legend.text = element_text(size= 12)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance\n")
3c - Phyla- and Proteobacterial class-level stacked barchart by DNA_or_cDNA~Month
tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")
ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")
4 - Top SVs and their sequences
Top20SVs_names = names(sort(taxa_sums(ps.final.r), TRUE)[1:20])
Top20SVs_names_df=as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6])
rownames(Top20SVs_names_df) = c()
Top20SVs_names_df## Phylum Class Order Family
## 1 Proteobacteria Alphaproteobacteria Rhizobiales Beijerinckiaceae
## 2 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 3 Proteobacteria Alphaproteobacteria Rhizobiales Labraceae
## 4 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 5 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 6 Proteobacteria Alphaproteobacteria Rhizobiales Xanthobacteraceae
## 7 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 8 Proteobacteria Alphaproteobacteria Rhizobiales Beijerinckiaceae
## 9 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 10 Proteobacteria Alphaproteobacteria Acetobacterales Acetobacteraceae
## 11 Actinobacteria Actinobacteria Micrococcales Microbacteriaceae
## 12 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 13 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 14 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 15 Actinobacteria Actinobacteria Micrococcales Micrococcaceae
## 16 Actinobacteria Actinobacteria Corynebacteriales Nocardiaceae
## 17 Actinobacteria Actinobacteria Frankiales <NA>
## 18 Actinobacteria Actinobacteria Micrococcales Cellulomonadaceae
## 19 Actinobacteria Actinobacteria Frankiales <NA>
## 20 Actinobacteria Actinobacteria Frankiales <NA>
## Genus
## 1 Bosea
## 2 Cupriavidus
## 3 Labrys
## 4 Sphingomonas
## 5 Burkholderia-Caballeronia-Paraburkholderia
## 6 Bradyrhizobium
## 7 Sphingomonas
## 8 Methylobacterium
## 9 Sphingomonas
## 10 <NA>
## 11 Salinibacterium
## 12 Massilia
## 13 Sphingomonas
## 14 Delftia
## 15 Arthrobacter
## 16 <NA>
## 17 <NA>
## 18 <NA>
## 19 <NA>
## 20 <NA>rownames(as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6]))## [1] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTGTCCGGGAAGATAATGACTGTACCGGAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGACTCTTAAGTCGGGGGTGAAAGCCCAGGGCTCAACCCTGGAATTGCCTTCGATACTGGGAGTCTTGAGTTCGGAAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGTCCGAAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [2] "TGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCGCTGGTTAATACCTGGCGTGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGACAGGCGTGAAATCCCCGGGCTTAACCTGGGAATTGCGCTTGTGACTGCAAGGCTAGAGTGCGTCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGACGTGACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [3] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGACGGCCTTAGGGTTGTAAAGCTCTTTTAACAGGGACGATAATGACGGTACCTGTAGAATAAGCCCCGGCAAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGAGGCTCAACCTCAGAACTGCCTTCGATACTGGCGATCTTGAGTTCGGAAGAGGTTGGTGGAACAGCTAGTGTAGAGGTGAAATTCGTAGATATTAGCTAGAACACCAGTGGCGAAGGCGGCCAACTGGTCCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [4] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCTGGTGCTCAACACCAGAACTGCCTTTTAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [5] "TGGGGAATTTTGGACAATGGGGGAAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAAACCTTTGCACTAATACTGTGAGGGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTCGTTAAGACAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGTGACTGGCGAGCTAGAGTATGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGCCAATACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [6] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTGTGCGGGAAGATAATGACGGTACCGCAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGTGCGTAGGCGGGTCTTTAAGTCAGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGAAGATCTTGAGTTCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGATACTGACGCTGAGGCACGAAAGCGTGGGGAGCAAACA"
## [7] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [8] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTATCCGGGACGATAATGACGGTACCGGAGGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGCGTTTTAAGTCGGGGGTGAAAGCCTGTGGCTCAACCACAGAATGGCCTTCGATACTGGGACGCTTGAGTATGGTAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGACCATCACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [9] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTACCCGGGAAGATAATGACTGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTCAGAGGTGAAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAGACTGCATCGCTTGAATCCAGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [10] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTCGGCGGGGACGATGATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATGACTGGGCGTAAAGGGCGCGTAGGCGGTTGCATAAGTTAGATGTGAAATTCCCGGGCTTAACCTGGGGACTGCATTTGATACTATGTGGCTTGAGTATGGAAGAGGGTCGTGGAATTCCCAGTGTAGAGGTGAAATTCGTAGATATTGGGAAGAACACCGGTGGCGAAGGCGGCGACCTGGTCCATAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [11] "TGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCAACGCCGCGTGAGGGACGACGGCCTTCGGGTTGTAAACCTCTTTTAGTAGGGAAGAAGCGAAAGTGACGGTACCTGCAGAAAAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTGGGGGCTCAACCCCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAATGGCGAAGGCAGATCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [12] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCAGGGAAGAAACGGTGTGGGCTAATATCCTGCACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCTGTCGTGAAATCCCCGGGCTCAACCTGGGAATTGCGATGGAGACTGCAAGGCTAGAATCTGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGTCAAGATTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [13] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTTAGAGGTGAAAGCCTGGAGCTCAACTCCAGAATTGCCTTTAAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [14] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGCAGGATGAAGGCCTTCGGGTTGTAAACTGCTTTTGTACGGAACGAAAAAGCTCCTTCTAATACAGGGGGCCCATGACGGTACCGTAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTATGTAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTTGTGACTGCATGGCTAGAGTACGGTAGAGGGGGATGGAATTCCGCGTGTAGCAGTGAAATGCGTAGATATGCGGAGGAACACCGATGGCGAAGGCAATCCCCTGGACCTGTACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [15] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGTAGGGAACAAGGCCACACGTGGTGTGGTTGAGGGTACTTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCCGTGAAAGTCCGGGGCTCAACTCCGGATCTGCGGTGGGTACGGGCAGACTAGAGTGATGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCATTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [16] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAAGCGAGAGTGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTGTCCGGAATTACTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCACGTCGGCTGTGAAAACCCATCGCTCAACGGTGGGCTTGCAGTCGATACGGGCTGACTTGAGTACTGCAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"
## [17] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"
## [18] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGGAGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGAAGAAGTGGTCGGGTTCTCTCGGTCGTGACGGTACCTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTCGAGGCTCAACCTCGGGCTTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCCGCAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [19] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCAAACA"
## [20] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAACACAGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTCTGTCGCGTCGGCTGTGAAAACTCGGGGCTCAACCCCGAGCCTGCAGTCGATACGGGCAGACTAGAGTGCTGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"5a - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (no filter)
5b - Genus-level stacked barchart by DNA_or_cDNA~Snowpack_Layer (filter)
Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.
ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)
ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE),
by = "Genus"]
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"
ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]
ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
levels = c("Acidiphilium","Arthrobacter",
"Bacillus",
"Bosea",
"Burkholderia-Caballeronia-Paraburkholderia",
"Cryobacterium","Cupriavidus",
"Delftia","Hymenobacter",
"Labrys","Leptothrix",
"Massilia","Methylobacterium",
"Noviherbaspirillum","Oryzihumus",
"Polymorphobacter",
"Salinibacterium",
"Sphingomonas","Other (<1%)"))tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
"#4477AA", "#77AADD", "#117777", "#44AAAA",
"#77CCCC", "#777711", "#AAAA44", "#DDDD77",
"#774411", "#AA7744", "#DDAA77", "#771122",
"#AA4455", "#DD7788", "azure3")
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Snowpack_Layer,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")
5c - Genus-level stacked barchart by DNA_or_cDNA~Fox_Aspect (filter)
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")
Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.
I don’t like this, but plotting all samples now, without grouping.
# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
# levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Sample_Name, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")
6 - Unconstrained ordination
nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
"June" = "#0704A6",
"Early July" = "#727A0B",
"Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
"Mid" = "#C6156D",
"Basal" = "#890909",
"Sup. Ice" = "#091D89",
"Glacial Ice" = "#580303")
ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")
ps.final.r.dna.nmds <- ordinate(ps.final.r.dna,
"NMDS",
distance="jsd",
maxit = 1e3)## Run 0 stress 0.2107623
## Run 1 stress 0.2116726
## Run 2 stress 0.2123577
## Run 3 stress 0.2133448
## Run 4 stress 0.2030457
## ... New best solution
## ... Procrustes: rmse 0.1022482 max resid 0.3790585
## Run 5 stress 0.2161223
## Run 6 stress 0.2165051
## Run 7 stress 0.2139681
## Run 8 stress 0.2087422
## Run 9 stress 0.1999096
## ... New best solution
## ... Procrustes: rmse 0.05625804 max resid 0.2574247
## Run 10 stress 0.2102544
## Run 11 stress 0.2150366
## Run 12 stress 0.2026495
## Run 13 stress 0.2154225
## Run 14 stress 0.2213739
## Run 15 stress 0.2098726
## Run 16 stress 0.2163894
## Run 17 stress 0.2031956
## Run 18 stress 0.2015388
## Run 19 stress 0.2044157
## Run 20 stress 0.212047
## *** No convergence -- monoMDS stopping criteria:
## 20: stress ratio > sratmaxps.final.r.cdna.nmds <- ordinate(ps.final.r.cdna,
"NMDS",
distance="jsd",
maxit = 1e3)## Run 0 stress 0.1654426
## Run 1 stress 0.1685203
## Run 2 stress 0.1630985
## ... New best solution
## ... Procrustes: rmse 0.02948222 max resid 0.1409178
## Run 3 stress 0.1644532
## Run 4 stress 0.1734149
## Run 5 stress 0.1665583
## Run 6 stress 0.1680772
## Run 7 stress 0.1626706
## ... New best solution
## ... Procrustes: rmse 0.120582 max resid 0.4752591
## Run 8 stress 0.1623674
## ... New best solution
## ... Procrustes: rmse 0.09400929 max resid 0.3635666
## Run 9 stress 0.166825
## Run 10 stress 0.1621902
## ... New best solution
## ... Procrustes: rmse 0.01076304 max resid 0.0346839
## Run 11 stress 0.161498
## ... New best solution
## ... Procrustes: rmse 0.06954823 max resid 0.2343931
## Run 12 stress 0.1714477
## Run 13 stress 0.1642582
## Run 14 stress 0.1663909
## Run 15 stress 0.1602385
## ... New best solution
## ... Procrustes: rmse 0.02810959 max resid 0.1054134
## Run 16 stress 0.1687445
## Run 17 stress 0.1691279
## Run 18 stress 0.1670851
## Run 19 stress 0.1690936
## Run 20 stress 0.1652288
## *** No convergence -- monoMDS stopping criteria:
## 17: stress ratio > sratmax
## 3: scale factor of the gradient < sfgrminps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.5, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.15,
linetype = 2,
show.legend = FALSE) +
ggtitle("DNA nMDS") +
guides(colour = guide_legend(title="Month"))
ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.7, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.2,
linetype = 2,
show.legend = FALSE)+
ggtitle("cDNA nMDS") +
guides(colour = guide_legend(title="Month"))plot_grid(ps.final.r.dna.nmds.plot + theme_bw(),
ps.final.r.cdna.nmds.plot + theme_bw())
7 - Alpha-diversity metrics.
month_comparisons = list( c("April", "June"), c("April", "Early July"),
c("April", "Late July"),
c("June", "Early July"), c("June", "Late July"),
c("Early July", "Late July"))
site_comparisons = list( c("NW_AWS", "N_NE"),
c("NW_AWS", "SWS1SE"),
c("N_NE", "SWS1SE"))
slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"),
c("TOP", "GL ICE"),
c("MID", "BASAL"), c("MID", "SUP ICE"),
c("MID", "GL ICE"),
c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
c("SUP ICE", "GL ICE"))All measures for Sample_ID, no stats, facet by DNA/cDNA
plot_richness(ps.final,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures by Sample_ID
estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>%
kable(digits = 4) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover")) %>% row_spec(0,bold=TRUE)| Observed | Shannon | Simpson | |
|---|---|---|---|
| A03.cDNA.1 | 286 | 0.7114 | 0.2103 |
| A03.DNA.1 | 138 | 3.2774 | 0.9407 |
| A04.cDNA.2 | 163 | 2.0361 | 0.7852 |
| A04.DNA.2 | 75 | 1.6409 | 0.6083 |
| A05.cDNA.3 | 515 | 2.2067 | 0.7048 |
| A05.DNA.3 | 1390 | 5.1230 | 0.9813 |
| A06.cDNA.4 | 125 | 1.3510 | 0.6058 |
| A06.DNA.4 | 210 | 3.9337 | 0.9650 |
| A07.cDNA.5 | 1635 | 4.6591 | 0.9488 |
| A07.DNA.5 | 609 | 3.3602 | 0.8726 |
| A08.DNA.6 | 563 | 4.0141 | 0.8786 |
| A09.cDNA.7 | 1875 | 5.5902 | 0.9807 |
| A10.cDNA.8 | 842 | 4.7339 | 0.9520 |
| A10.DNA.8 | 179 | 3.3916 | 0.9493 |
| A11.cDNA.9 | 3 | 0.2518 | 0.1121 |
| A11.DNA.9 | 1819 | 5.3121 | 0.9764 |
| A12.cDNA.10 | 54 | 1.4929 | 0.6133 |
| A12.DNA.10 | 2719 | 5.7179 | 0.9837 |
| B01.cDNA.11 | 453 | 4.0612 | 0.9622 |
| B01.DNA.11 | 187 | 3.5277 | 0.9546 |
| B02.cDNA.12 | 733 | 4.5851 | 0.9683 |
| B02.DNA.12 | 310 | 2.1151 | 0.7154 |
| B03.cDNA.13 | 38 | 1.1241 | 0.4091 |
| B03.DNA.13 | 457 | 4.2859 | 0.9559 |
| B04.cDNA.14 | 67 | 1.1867 | 0.4139 |
| B05.cDNA.15 | 1307 | 5.7881 | 0.9881 |
| B06.cDNA.16 | 134 | 1.4982 | 0.6595 |
| B06.DNA.2 | 156 | 0.6586 | 0.2271 |
| B07.DNA.3 | 1302 | 5.3971 | 0.9910 |
| B08.cDNA.18 | 21 | 0.9250 | 0.5154 |
| B08.DNA.4 | 986 | 5.2137 | 0.9805 |
| B09.DNA.5 | 105 | 2.9238 | 0.9150 |
| B10.cDNA.21 | 400 | 4.5700 | 0.9771 |
| B10.DNA.6 | 377 | 4.1547 | 0.9598 |
| B12.cDNA.23 | 29 | 0.9625 | 0.3511 |
| B12.DNA.8 | 1998 | 5.3382 | 0.9817 |
| C01.DNA.9 | 165 | 3.3370 | 0.9306 |
| C02.cDNA.26 | 907 | 4.1753 | 0.9419 |
| C03.cDNA.27 | 537 | 3.3393 | 0.8840 |
| C03.DNA.B2.1 | 677 | 5.1429 | 0.9887 |
| C04.cDNA.28 | 307 | 3.3747 | 0.8749 |
| C04.DNA.2 | 715 | 4.2502 | 0.9430 |
| C05.DNA.3 | 1739 | 5.2347 | 0.9838 |
| C06.DNA.4 | 236 | 3.4402 | 0.9006 |
| C07.DNA.5 | 2789 | 5.9749 | 0.9894 |
| C08.DNA.6 | 307 | 4.5001 | 0.9837 |
| C09.cDNA.33 | 36 | 1.3284 | 0.5802 |
| C09.DNA.7 | 148 | 3.8243 | 0.9690 |
| C10.DNA.8 | 383 | 3.2676 | 0.8801 |
| C11.cDNA.35 | 65 | 1.5747 | 0.6657 |
| C12.cDNA.36 | 19 | 0.4213 | 0.1462 |
| C12.DNA.10 | 396 | 2.2692 | 0.7492 |
| D01.DNA.B4.1 | 147 | 0.6873 | 0.2422 |
| D02.cDNA.38 | 128 | 2.1030 | 0.7780 |
| D03.cDNA.39 | 266 | 3.4007 | 0.9035 |
| D04.cDNA.40 | 79 | 1.5083 | 0.6114 |
| D05.cDNA.41 | 154 | 1.8281 | 0.7377 |
| D05.DNA.5 | 90 | 0.6634 | 0.2169 |
| D06.cDNA.42 | 841 | 3.6042 | 0.9024 |
| D06.DNA.6 | 2292 | 5.5649 | 0.9810 |
| D07.cDNA.43 | 1444 | 4.5026 | 0.9489 |
| D07.DNA.7 | 467 | 1.8746 | 0.5308 |
| D08.cDNA.44 | 120 | 2.6011 | 0.8428 |
| D08.DNA.8 | 1045 | 4.1510 | 0.9245 |
| D09.DNA.9 | 2986 | 5.9104 | 0.9895 |
All measures for Month, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Month, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Month, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
mini_month_comparisons = list(c("April", "June"),
c("April", "June"))
plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Fox_Aspect, no stats, facet by DNA/cDNA
#now without stats, only shapes
plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Fox_Aspect, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed
plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, no (“Early July”, “Late July”)
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month != "Early July" & Month != "Late July"),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, July only, “GL ICE” removed, “Gl ice” samples removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” also removed."
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"),
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Early July” only, “GL ICE” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
All measures for Fox_Aspect, stats included, facet by DNA/cDNA, “Late July” only, “GL ICE”, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
All measures for Snowpack_Layer, no stats, facet by DNA/cDNA
#now without stats, only shapes
plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Snowpack_Layer, stats included, facet by DNA/cDNA
plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA
plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)All measures for Snowpack_Layer in July only (“Early July”, “Late July”), stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed, “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_D” removed
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
Snowpack_Layer != "GL ICE" &
!(Sample_Name %in% gl_ice_spls) &
Sample_Name != c("30 July_SWS1SE_Slush_R",
"30 July_SWS1SE_Slush_D")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)All measures for Snowpack_Layer in “Early July”, stats included, facet by DNA/cDNA, “GL ICE” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)All measures for Snowpack_Layer in Late July, stats included, facet by DNA/cDNA, “GL ICE” removed, samples “30 July_SWS1SE_Slush_R”, “30 July_SWS1SE_Slush_R” removed, “Gl ice” samples removed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning in max(data$density): no non-missing arguments to max; returning -Inf## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)8 - Alpha diversity with rarefaction to the smallest number of reads in any sample
ps.final.rare = rarefy_even_depth(ps.final,
sample.size = min(sample_sums(ps.final)),
rngseed = TRUE, replace = TRUE,
trimOTUs = TRUE, verbose = TRUE)## `set.seed(TRUE)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(TRUE); .Random.seed` for the full vector## ...## 25951OTUs were removed because they are no longer
## present in any sample after random subsampling## ...plot_richness(ps.final.rare,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
plot_richness(ps.final.rare,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Month") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
#now without stats, only shapes
plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
#now without stats, only shapes
plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))
plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))